Universality in Fluid Domain Coarsening: The case of vapor-liquid transition 
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Domain growth during the kinetics of phase separation is studied following vapor-liquid transition 
in a single component Lennard- Jones fluid. Results are analyzed after appropriately mapping the 
continuum snapshots obtained from extensive molecular dynamics simulations to a simple cubic 
lattice. For near critical quench interconnected domain morphology is observed. A brief period of 
slow diffusive growth is followed by a linear viscous hydrodynamic growth that lasts for an extended 
period of time. This result is in contradiction with earlier inclusive reports of late time growth 
exponent 1/2 that questions the uniqueness of the non-equilibrium universality for liquid- liquid and 
vapor-liquid transitions. 

PACS numbers: 64.60.Ht, 64.70.Ja 



I. INTRODUCTION 

Understanding non-equilibrium evolution during the 
phase separation in a system is of fundamental importance 
[l|-[4[ and of much research interest, both theoretically and 
experimentally [5l-[22j. Upon quenching from a high tem- 
perature homogeneous state to a temperature below the 
critical point, the system becomes unstable to fluctuations 
and starts phase separating with the formation and usu- 
ally nonlinear growth of particle rich and particle depleted 
domains. Such coarsening is a scaling phenomenon [23], 
e.g., the shape functions characterizing the morphology 
obey the scaling relations [3-0? Ell Hil 



C(r,t) = C(r/£(t)), 
S(k,t) = £(t) d S(ke(t)), 
P(£ d ,t) = t{t)- l P[l d /£{t)}, 



(1) 
(2) 
(3) 



where C(r, i), S(k, t) and P(£d : i) are respectively the two- 
point equal-time correlation function, structure factor and 
domain size distribution function. In Eqs. (p3|3]) C(x), 
S(y) and P(z) are the master functions independent of 
time (t)-dependent average domain-size £{t) that typically 
grows in a power law manner 



£(t) ~ t° 



(4) 



The growth exponent a, in Eq. (j4]), depends upon the 
transport mechanism driving the phase separation. For 
diffusive dynamics, associating the domain growth with 
the chemical potential gradient as 



d£(t) 
dt 



Ml 



7 



£{tf 



(5) 



7 being the inter- facial tension, one obtains a = 1/3. This 
is referred to as the Lifshitz-Slyozov [24[ (LS) law. While 
LS law is the only expected behavior for phase separating 
solid mixtures, one expects much faster growth, at large 
length scales, for fluids and polymers, where hydrodynam- 
ics is important. Compared to the kinetics in vapor-liquid 
phase separation, domain coarsening in binary fluids re- 
ceived much more attention [lE [13, [H, [H, SH-EH 



where following consensus for the behavior of late stage 
growth [29l432| has been reached in d = 3. Considering a 
balance between the surface energy density ( / y/£(t)) and 
the viscous stress (6irrjve/£(t), V£ being the interface veloc- 
ity and 77 the shear viscosity) for an interconnected domain 
structure, one can write 



d£(t) 
dt 



v £ -. 



(6) 



Solution of Eq. (0 predicts a linear growth (a = 1), a 
picture that holds for low Reynolds number. However, for 
£{t) ^> £i n (= ^ 2 /(p7), p being the density), the inertial 
length, the surface energy density is balanced by kinetic 
energy density pvf, so that 



d£(t) 
dt 



1 



£{ty/ 2 ' 



(7) 



giving a = 2/3. Note that a = 1 is referred to as the 
viscous hydrodynamic growth and a = 2/3 as inertial 
hydrodynamic growth. While crossover from diffusive to 
viscous hydrodynamic regime was observed in molecular 
dynamics (MD) simulations [HI, [HI and experiments [33l — 
[33, both viscous and inertial growths were observed in 
lattice-Boltzmann simulations [3|- 

While non-equilibrium universality for vapor-liquid 
phase separation is expected to be the same as the liquid- 
liquid, rare inclusive MD simulations that exist |15l. I36i - l39j 
for the vapor-liquid transition report an exponent a = 1/2 
for late time domain coarsening. Thus, even though our 
recent focus has turned to systems with realistic interac- 
tions and physical boundary conditions (EI-HH. our basic 
knowledge and understanding of segregation kinetics in 
bulk fluids still remains a challenge. In this letter, we 
present results from large scale MD simulations for kinet- 
ics of vapor-liquid phase separation in a single-component 
Lennard- Jones (LJ) fluid. This work reports first observa- 
tion of viscous hydrodynamic growth via MD simulation, 
thus computationally confirming that the non-equilibrium 
universality class of vapor-liquid and liquid-liquid phase 
separation is indeed the same as in the case of equilibrium 
critical phenomena for both static and dynamic proper- 
ties. 
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FIG. 1: (a) Left frames show the 3 — d snapshots during evolu- 
tion of the vapor-liquid phase separation in a simple Lennard- 
Jones system, at two different times, as indicated. The total 
density was set to 0.3 that gives interconnected domain mor- 
phology,as seen. The simulation was performed after quench- 
ing a homogeneous system prepared at a high temperature to a 
temperature T = 0.7. The right panel corresponds to the same 
snapshots, but after removing the noise (see text for details) 
to obtain a pure domain morphology, (b) Cross-sectional view 
of the same system in xy plane at two other times. Here also 
the left frames show the original snapshots and the right ones 
show the same after removing noise. 

II. MODEL AND METHOD 

We consider a model where particles interact via 
Vim) = U{r l3 ) - U(r c ) - (r« - r c ) (^^) ,(8) 

V l 3 / nj—r c 

with 



" / x 12 








\ r ij) 





being the standard LJ potential. In Eqs. J5J) and (|9j) 
is the scalar distance between the zth and jth particles, a 
is the particle diameter and e is the interaction strength. 
The cut-off and shifting of the potential at r c was used to 
faclitate faster computation. Further, introduction of the 
third term [So| on the right hand side ensures that both 
the effective potential and force are continuous at r c which 
was set to 2.5<r. All particles were assigned equal mass m. 
For the sake of convenience, below we set m, a, e and ks 
to unity. 

The phase-diagram of the model without the force- 
correction term is known [41], the critical temperature 
(T c ) and critical density (p c ) being 1.085 ± 0.005 and 
0.317±0.006 respectively. It is understood, of course, that 
the introduction of force-correction will somewhat change 
the value of T c , which, from our experience [42-44] with 
a similar model that was used to study phase-separation 
in binary mixture, could be approximately 10% lower. In 
this work we have chosen the temperature low enough so 
that even without the complete knowledge of the critical 
parameters, we are safely below the co-existence curve, as 
is also clear from the nice phase-separating domain pat- 
tern shown in Fig. [IJ The interconnected nature of the 
morphology also ensures that our choice of overall density 
(p) is close to the critical value. In a future correspon- 
dence, of course, we will undertake the understanding of 
the phase behavior of this model in conjunction with the 
presentation of detailed results on the phase-separation 
kinetics at various temperatures and densities. 

A total 265421 particles were confined in a periodic 
3-dimensional box of linear dimension L = 96. Thus we 
worked with a density p = 0.3. The MD runs were per- 
formed using the Verlet velocity algorithm (SO, EH with an 
integration time step of At = 0.005r, r = (ma 2 /e) 1 / 2 = 1. 
The temperature T was controlled by using a Nose- Hoover 
thermostat (NHT) [45] , which is known to preserve hydro- 
dynamics well. At t = homogeneous initial configura- 
tions, prepared by equilibrating the system at T = 5, were 
quenched to T < T c . All the results were obtained by av- 
eraging over 5 independent runs at a quench temperature 
T = 0.7. 

The left frames of Fig. []Ja) are typical 3-dimensional 
evolution pictures from two different times where the black 
dots represent a particle. They demonstrate nice inter- 
connected structure of segregating domains growing with 
time. The average density of liquid and vapor domains 
are respectively 0.72 and 0.02. A direct method to char- 
acterize the structure is to calculate the radial distribu- 
tion function, g(r), from these continuum configurations. 
However, we are interested in looking at C(r, i), its Fourier 
transform S(k, t) and P{1^ t) which have been more stan- 
dard tools in the study of phase ordering dynamics. To 
facilitate the calculation of latter quantities, we map the 
continuum snapshots to a lattice one in such a way that a 
pure domain structure is obtained where all lattice sites in- 
side a liquid domain have occupancy one while sites inside 
the vapor domain have occupa ncy zero. An appropriate 
renormalization procedure [HI, [l6|, [l9| is applied for this 
noise removal exercise. We further map this system into 
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spin- 1/2 Ising system where an occupied lattice site has 
spin value Si = 1 and for a vacant site Si = — 1 . The right 
frames in Fig. [TJa) show such mapped configurations. 
In Fig. DIb) we demonstrate similar exercise by showing 
cross-sectional views at two other times. Which may be 
difficult to appreciate from the 3 — d snapshots, these 2 — d 
pictures clearly show that no significant structural change 
has occurred during this transformation that can hamper 
the analysis. Also during the entire duration of study, we 
find that the composition of up and down spins remains 
conserved within a tolerable limit of 5% fluctuations. 

From the mapped configurations, the correlation func- 
tion C(r) is calculated as 



C(r,t) = (SiSj) - (Stf 



(10) 



where the angular brackets stand for statistical averaging. 
The average domain size £(t) from the decay of C(r, t) was 
estimated as 



C(r = £(t),t) = h, 



(11) 



where h was set to the first zero of C(r, t). Other methods 
employed for the calculation of £(t) are 



£(t) 



1 



Jdk kS(k) 



and 



e(t) = J de d e d p(e d ,t). 



(12) 



(13) 



In Eqs. (fT2j) and ([T3]h it is assumed that S(k,t) and 
P(£d,t) are appropriately normalized. In Eq. (fl3j) . £d 
is obtained from the separation between two successive 
domain interfaces in x—, y— or z— directions. 



III. RESULTS 

In Fig. [2ja), we present the scaling plot of C(r, t) vs. 
r/£(t) for which four different times were chosen, as indi- 
cated on the figure. While collapse of data from the ear- 
liest time with the other ones is not good, excellent data 
collapse with three later times is indicative that a scaling 
regime is reached. In Fig. Efb) we show the scaling behav- 
ior of S(k, t) on a log-log plot. The parallel nature of the 
tail region to the continuous line confirms consistency with 
the Porod law (461,113. Fig. Efc) shows the scaling plot of 
P(£d, t) on a semi- log plot where the linear behavior of the 
tail region is consistent with an exponential decay [l6|, [20| . 
Note that in all the cases, £(t) used was obtained from Eq. 
(fT3|) . The nice data collapse obtained in all the quantities 
by using £(t) from a single measure speaks about the con- 
sistency of different methods. Henceforth all the results 
for £(t) will be presented after calculating via Eq. ([T3|) . 

Coming to the central objective of the paper, viz., quan- 
tifying the dynamics, in Fig. [3] we present plot of £{t) vs. 
f on a log scale. There various lines correspond to var- 
ious growth laws in domain coarsening problem. While 
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FIG. 2: Scaling plot of (a) C(r,t) vs r/£(t) (b) £(t)~ d S(k,t) 
vs k£(t) and (c) £(t)P(£d,t) vs £d/£(t) where data from four 
different times, as indicated, have been collapsed. All quan- 
tities were calculated from noise-free domain structures. The 
continuous line in (b) corresponds to the Porod law while the 
one in (c) represents an exponential decay of the tail. 



the early time data is more consistent with a diffusive 
growth (a = 1/3), there appears a gradual crossover to 
a faster growth. The intermediate time data, of course, 
look consistent with t 1 ^ 2 as reported in previous works 
[HI, l36l - [39| and at a later time there seems to be a further 
crossover. However, due to large offset at the crossover 
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FIG. 3: Plot of £{t) vs t on a log-scale. All data were obtained 
after removing the noise as described in the text. Here different 
lines correspond to growth laws £{t) ~ t 1/3 , t 1/2 , t 2/3 and t, 
as indicated. 



region(s), which will be clearer from further analysis and 
discussions, we warn the reader not to take this conclusion 
seriously. For similar reasons, an apparently t 1 / 3 looking 
behavior during early time should also not be taken seri- 
ously. In fact, due to this and other practical difficulties, 
it is advisable not to reach a firm conclusion from looking 
at the behavior on a log-scale (which, of course, is a stan- 
dard practice followed in the literature), unless one has 
results over several decades in time and length. 
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FIG. 4: Plot of £(t) vs t on a linear scale. All data were 
obtained after removing the noise as described in the text. 
The solid line corresponds to the expected linear growth in 
the viscous hydrodynamic regime. Inset: Plot of the growth 
exponent a, as a function of 1/to, obtained by fitting the data 
in (a) to the form (fT4|) in the range [to, 350]. 



Fig. H] shows a plot of £(t) vs t on a linear scale where 
the data starting from t ~ 100 look quite linear all the 
way till the end. The deviation from this linear behavior at 
early time (t < 100) could be attributed to slower diffusive 
growth, as discussed above. A fitting to the form, 



£{t) = A + Bt a , 



(14) 



however, in the range t e [0, 100] gives a = 0.66. Since 
this diffusive regime is short lived and also accompanied 
by gradual crossover to a faster growth regime, it is in 
fact not practical to search for growth exponent a — 1/3 
in this region. On the other hand, a similar fitting to the 
range t e [100,350] gives a ~ 0.9, quite consistent with 
the predictions for viscous hydrodynamic growth. Note 
that for t > 350 one hits the finite-size effects for the 
system size used in this work. In the inset of Fig. |4] 
we present plot of the exponent a as a function of 1/to 
obtained by fitting the data to the form ([14]) in the range 
[to, 350]. This exercise conclusively says that the exponent 
at late time certainly is larger than 1/2 which was reported 
earlier. Even though in the main frame of Fig. 3] one sees 
a linear behavior for an extended period of time, for the 
sake of completeness and to further strengthen our claim 
we take the following route which, in addition to being a 
nice exercise, could be useful for other complex situations. 

In Fig. [5] we present a plot of instantaneous exponent 
oti calculated as [48[ 



d[\n£(t)} 
d[\nt] 



(15) 



vs 1 /£(£). Due to the off-set (^o) at t = 0, as discussed in 
the context of Fig. [3j oti is less than unity for the whole 
range. However it seems to be approaching unity in a 
non-linear fashion. The dashed line there with an arrow 
at the end serves as a guide to the eyes. The lower part 
in this figure represents the corresponding result obtained 
from MD simulation by using Andersen thermostat (AT) 
[45| . Due to the stochastic nature of this thermalization 
algorithm, it is expected that a diffusive growth will be 
seen for all time and lengthscales and indeed is seen. There 
the dashed line has the form [12| 



£(t) 



(16) 



which says that the observation of a < 1/3 for finite do- 
main lengths is a mathematical artifact coming due to 
non-zero value of £o and does not necessarily mean that the 
actual exponent is not 1/3 at early time. Coming to the 
point of hydrodynamic preserving capability of NHT, of 
course, much better thermostats are available these days. 
However, from the difference seen for the results coming 
from NHT and AT, it is clear that NHT is preserving it 
rather well. Thus the validity of the methodology adopted 
in this work is justified. 

Note that in Eq. (fTBj) £q corresponds to the character- 
istic length scale at the beginning of a scaling regime and 
not necessarily the length at t = 0. In situations where 
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FIG. 5: Plot of instantaneous exponent on vs l/t(t). The 
dashed line is a guide to the eyes, while the dotted line cor- 
responds to Eq. ([T6|) . Inset: Plot of on vs l/i'(t') for two 
different values of t c , as indicated. There the horizontal line 
represents there the linear hydrodynamic growth. 



one does not expect multiple scaling regimes (as is the case 
with AT) £o can be nicely estimated from a data-collapse 
experiment in finite-size scaling analysis [l2|, |49| . But situ- 
ation is more complex with NHT due to a crossover to hy- 
drodynamic regime and non- availability of data for differ- 
ent system sizes. However, considering that the crossover, 
at time t c , under discussion is to a linear regime, confir- 
mation about it could be obtained by using the following 
simpler method. To start with, we assume that the growth 
obeys a power law behavior as a function of the shifted 
time t' = t — t c as 



£'(t')=£(t)-£(t c )=At' a 
and calculate the instantaneous exponent 

d[\n£'(t')} 



d[lnt'] 



(17) 



(18) 



Eq. ([TT]) is invariant under an arbitrary choice of t c if we 
are in a linear growth regime and one should obtain a con- 
stant value oti — 1 for all choices of t c in the post crossover 
regime. In the inset of Fig. [5j we present oti for two dif- 
ferent choices of t c . For t c > 110, oscillates around 
the mean value 1 which is consistent with the viscous hy- 
drodynamic regime. This oscillation (which could as well 



be appreciated from main part of Fig. [5j) is also observed 
in other studies [HI, HH • In a finite system, at later time 
only a few domains of comparable size exist which are sep- 
arated from each other by large distances. Thus, as time 
increases, it takes longer time for them to merge which 
causes the exponent to come down. Finally merging of 
two huge domains after a long interval suddenly enhances 
the value of the exponent. 

IV. CONCLUSION 



We have presented results for the kinetics of vapor- 
liquid phase-separation from the molecular dynamics 
(MD) simulation of a simple Lennard- Jones fluid us- 
ing both Andersen (AT) and Nose-Hoover (NHT) ther- 
mostats. It is observed that NHT is reasonably useful for 
studying hydrodynamic effects in the fluid phase separa- 
tion. A brief period of diffusive coarsening was followed 
by a linear viscous hydrodynamic growth. Our results 
are in contradiction with few previous MD studies which 
reported an exponent much less than unity, however, is 
consistent with the results of binary fluid phase separa- 
tion. One requires much larger system size to observe 
the inertial hydrodynamic growth. This we leave out for 
a future exercise, in addition to the study of coarsening 
of liquid droplets for off-critical quench. For the latter 
problem one expects single asymptotic exponent 1/3 from 
droplet diffusion-coagulation mechanism [23, 50]. Also, as 
discussed, MD with AT and a Monte Carlo simulation of 
the same system will provide a single diffusive 1/3 growth 
for all compositions. For off-critical composition the lat- 
ter exercise should be able to provide information on the 
curvature dependent inter-facial tension which is expected 
to bring in early-time correction in the LS law. This will 
certainly be interesting to compare with the corresponding 
results obtained from equilibrium studies [51| . 
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